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I. 


INTRODUCTION 


In  the  last  few  years,  a concerted  effort  has  been  made  to  analyze 
submerged  structures  under  explosive  loadings,  by  utilizing  the  Doubly 
Asymptotic  Approximation  (DAA) , Ref.  [1],  to  account  for  the  structure-fluid 
interaction.  The  current  implementation  of  the  DAA.by  Weidlinger  Associates, 
Ref.  [2]  to  [5],  Employs  a modal  analysis  for  the  structure  and  a set  of 
orthogonal  functions  (surface  expansion  functions)  for  the  fluid,  the  latter 
being  used  to  improve  the  conditioning  of  the  governing  response  equations. 

The  DAA  yields  exact  results  for  high  and  low  frequencies  and  produces  a 
smooth  transition  between  these  limits.  Since  a more  accurate  fit  at 
intermediate  frequencies  may  prove  to  be  important  in  some  problems  of 
interest,  the  DAA  may  prove  to  be  inadequate  in  some  cases.  Thus,  work  was 
done  on  developing  improvements  to  the  DAA.  One  such  approach,  the  Inertial- 
Damping  Collocation  Approximation  (IDCA),  was  developed  by 

Weidlinger  Associates  several  years  ago  and  was  reported  on  in  Ref.  [2].  The 
IDCA  gives  exact  results  for  low,  high,  and  selected  intermediate 
frequencies,  and  as  such  represents  a potential  improvement  over  the  DAA. 

Although  the  IDCA  was  described  in  Ref.  [2],  problems  arising  from  the 
fitting  procedures  used  to  match  exact  steady-state  results  were  not 
discussed  in  detail.  With  some  current  interest  in  the  development  of 
improved  uncoupling  techniques,  this  note  has  been  prepared  to  document  the 
current  status  of  the  IDCA  and  to  demonstrate  the  new  improved  fitting 
procedures  developed  for  the  application  of  this  method. 
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II.  DESCRIPTION  OF  THE  INERTIAL-DAMRINC  COLLOCATION  APPROXIMATION 


For  a submerged  structure  subjected  to  an  underwater  explosion,  it  is 
possible  to  approximate  the  structure-fluid  interaction  by  expressing  the 
surface  piessures  in  terms  of  the  normal  surface  velocities  (Refs.  [1]  and 
[2]).  Such  a procedure  essentially  uncouples  the  fluid  from  the  structure 
by  subjecting  the  structure  in-vacuo  to  free-field  effects  plus  a structure- 
fluid  interaction  loading.  The  DAA,  for  example,  is  an  uncoupling  technique 
which  yields  correct  results  for  low  and  high  frequencies.  In  this  section, 
a new  uncoupling  scheme,  the  IDCA,  will  be  described  briefly.  This  method 
allows  matching  exact  solutions  for  low,  high,  and  selected  intermediate 
frequencies . 


Let  a time-harmonic  normal  displacement  be  applied  to  a fluid  surface 
(cavity)  having  the  shape  of  a structure  to  be  submerged  and  loaded  by  an 
underwater  explosion: 


w(s,  t)  = W ( s ) e 


where  s denotes  position  on  the  surface,  t represents  time,  and  i 
The  radiated  fluid  pressure  may  then  be  written  in  the  form 


‘ -sF^ 


p(s,  t)  = [ p ( s ) + i p(s) 


in  which  the  inertial  (real)  component,  p,  is  in  phase  with  the 


displacement,  and  the  damping  (imaginary)  component,  p,  is  in  phase  with  the 
velocity. 


For  arbitrary  excitation  of  a submerged  structure,  the  radiated 
pressure  (the  difference  between  the  total  and  incident  pressures  on  the 
surface  of  the  structure)  may  also  be  represented  by  a component  in  phase 


with  the  displacement  of  the  surface  of  the  structure  and  a component  in 


phase  with  the  normal  velocity  of  the  surface.  The  generalized  force 


associated  with  the  radiated  pressure  may  then  be  expressed  as 

9 = 9j  + 92  (3) 

where  and  have  the  characters  of  an  inertia  force  and  a damping  force, 
respectively,  underlining  denotes  a matrix,  and  the  caps  indicate  a 
formulation  of  the  fluid  response  based  on  surface  expansion  functions  rather 
than  structural  modes. 


When  based  on  surface  expansion  functions,  the  IDCA  may  be  written  as 


*) 


t t 


Uj  dt  - a (q  - U^.)  = y 


-1 


Q dtdt  + B Q + C Q 


(4) 


q - Ui  = D 


§2dt  “ - 1 §2 


(5) 


in  which  q is  the  array  of  generalized  coordinates  for  the  fluid  , 

Uj  is  the  array  of  expansion  coefficients  for  the  incident 

/V  A 

(free-field)  velocity,  y is  the  virtual  mass  matrix,  a is  the  diagonal 
radiation  damping  matrix,  and  a dot  indicates  differentiation  with  respect  to 
time.  The  elements  of  matrices  B,  C,  D are  evaluated  by  matching  exact 
steady-state  solutions  at  selected  frequencies. 


For  low  frequencies,  Eq . (4)  reduces  to  the  Virtual  Mass  Approximation 
(Ref.  [6]),  and  for  high  frequencies,  Eq . (5)  reduces  to  the  Plane  Wave 
Approximation  (Ref.  [7]).  The  matrix  a in  Eq . (4)  is  used  to  ensure  the 

*j  " ~ ~ 

If  a submerged  body  is  subjected  to  an  incident  pressure  wave  with  zero 
rise  time,  a time-integrated  form  of  Eqs.  (4)  and  (5)  should  be  used, 
as  shown  in  Ref.  [2]. 
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correct  asymptotic  behavior  of  the  inertial  part  of  the  IDCA  for  high 
frequencies  (Curved  Wave  Approximation,  Ref.  [8]).  Details  of  the  procedures 
used  to  obtain  the  elements  of  a,  B,  C,  D may  be  found  in  Ref.  [2]. 
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III.  DIFFICULTIES  ENCOUNTERED 


In  the  IDCA,  fitting  is  employed  for  the  relations  between  normal 
surface  displacements  and  the  resulting  generalized  forces  considered  as 
functions  of  the  frequency.  The  concept  of  fitting  was  developed  after  a 
study  of  simple  cases  (infinitely  long  cylindrical  and  spherical  shells) 
where  one  displacement  mode  and  one  force  are  fitted  at  a time.  In  this 
situation  the  fitted  curves  are  necessarily  smooth,  and  in  general  acceptable. 
In  the  general  case,  a considerable  number  of  surface  expansion  coordinates 
and  the  corresponding  forces  are  fitted  simultaneously,  and  the  success  of 
this  fitting  process  was  a major  matter  of  concern  for  the  investigators. 

Discussing  as  a typical  example  the  case  N = 1 (one  circumferential 
wave)  for  the  structure  treated  in  Ref.  f 3 ] (ring-stiffened  circular 

cylindrical  shell  closed  by  flat  end  plates),  it  was  found  that  the  fitting 

*) 

process  properly  reproduced  the  desired  values  of  the  fluid  resistance  at 
the  frequencies  selected  for  fitting.  However,  the  curves  were  not  as  smooth 
as  hoped.  While  the  oscillations  in  most  cases  are  minor,  there  are  a number 
of  cases  with  violent  oscillations  in  narrow  frequency  ranges.  In  the  case 
under  discussion,  18  surface  expansion  functions  and  the  corresponding  18 
generalized  forces  were  used,  and  the  resulting  relations  plotted.  In  two 
cases  the  values  Y (accession  to  mass)  oscillated  so  severely  that  the 
values  became  negative  for  a narrow  frequency  range.  Figures  1 to  4 show  as 


In  this  note,  the  term  fluid  resistance  is  used  to  denote  the  array 

A A A A 

y = yd  + 1 Yt » which  YD  and  YT  have  real  elements  and  i = W-  1 . 

R L R I 

A 

The  array  YD  has  the  character  of  the  specific  acoustic  inertia,  while 

fl  Yj>  where  (1  represents  frequency,  has  the  character  of  the  specific 

acoustic  resistance.  The  exact  frequency  dependence  of  Y was  obtained 
by  modifying  an  existing  sound  radiation  code  (Ref.  [9]). 


r 
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A 

examples  plots  of  y for  surface  expansion  functions  1,  2,  5 and  12.  For 
R 

function  1,  the  rigid  body  displacement,  y is  in  general  well  approximated, 

R 

but  there  are  slight  oscillations  below  f - 0.5  kcps.  For  functions  2 and 
12,  Figs.  2 and  4 show  the  worst  situations  encountered;  they  give 
unrealistic  negative  values  for  y at  certain  frequencies.  An  intermediate 
situation  for  function  5 is  shown  in  Fig.  3.  There  are  strong  oscillations, 

A A 

but  no  negative  values  for  y . The  fitting  of  the  values  y (damping)  does 
not  lead  to  violent  oscillations  as  was  the  case  in  Figs.  2 to  4.  The 
result  for  the  rigid  body  motion,  function  1,  is  shown  in  Fig.  5.  However, 
inspection  of  the  results  for  the  remainder  of  the  surface  expansion 
functions  disclosed  two  cases,  functions  12  and  13,  shown  in  Figs.  6 and  7, 
where  the  damping  y becomes  negative  for  frequencies  below  f ~1.5  kcps. 

It  would  be  legitimate  to  reason  that  local  oscillations  in  the  curves 

for  y and  y should  not  cause  unrealistic  final  results,  because  the 
R L 

response  is  inherently  a composite  of  effects  originating  from  all 
frequencies.  However,  the  occurrence  of  regions  of  negative  damping  or 
negative  accession  to  mass  makes  the  results  of  the  analysis  tainted, 
although  only  4 of  the  total  of  18  relations  are  thus  affected.  Thought  was 
therefore  given  to  variations  of  the  scheme  which  avoided  the  described 


undesirable  situations. 
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IV.  IMPROVED  FITTING  PROCEDURE 


If  the  IDCA  in  terns  of  surface  expansion  functions,  as  described  in 

Ref.  [2],  is  applied  to  an  infinitely  long  cylindrical  shell  (or  to  a 

spherical  one)  of  uniform  thickness,  the  difficulties  described  in 

Section  III  could  not  occur,  because  only  one  surface  expansion  coordinate 

and  one  generalized  force  couple.  This  simple  situation  is  caused  by  the 

fact  that  the  matrices  y and  YT  are  diagonal.  This  suggests  that  one  could 

avoid  these  difficulties  if  one  could  select  surface  expansion  functions 

such  that  the  matrices  y and  y are  diagonal  at  all  frequencies.  This 

K 1 

ideal  solution  is  not  available  for  a general,  complex  structure. 


However,  one  can  start  by  selecting  surface  expansion  functions  i^.(s) 

A A 

which  make  the  matrices  y and  y diagonal  .for  ft  = 0 and  ft  °°, 
respectively,  where  ft  denotes  frequency.  This  is  always  possible,  and  can 
be  done  as  described  in  the  Appendix.  (Note  that  the  previously  proposed 
functions  tjjj(s)  already  lead  to  a diagonal  matrix  y^  for  ft  -»  °°.)  At 
intermediate  frequencies,  such  as  those  selected  for  fitting,  the  matrices 

A A A A 

y and  y would  not  be  diagonal.  However,  since  y and  y are  diagonal  for 
R I R X 

ft  = 0 and  ft  -»•  00 , respectively,  one  may  expect  that  the  matrices  at 
intermediate  frequencies  are  sufficiently  close  to  diagonal  to  permit 
ignoring,  as  an  approximation,  the  off-diagonal  terms.  Such  an  approach 
leads  to  an  IDCA-solution  without  problems  caused  by  negative  accession  to 
mass  or  negative  damping. 


In  order  to  evaluate  this  idea,  a unique,  but  of  course  truncated  set 

A 

of  surface  expansion  functions  was  obtained,  such  that  the  matrices  y and 

R 

Yj  were  diagonal  for  ft  = 0 and  ft  -*•  00 , respectively.  This  was  again  done 
for  the  case  studied  in  Section  III.  Using  these  new  surface  expansion 


functions,  the  matrices  y and  y were  computed  at  several  intermediate 

R i 


A 
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frequencies  to  verify  the  order  of  magnitude  of  the  off-diagonal  terms  to  be 
ignored . 


In  the  earlier  scheme  using  the  complete  matrix  of  the  fluid  resistance, 
including  the  off-diagonal  terms,  certain  fitting  frequencies  were  selected. 
Obviously,  the  importance  of  the  off-diagonal  terms  to  be  dropped  in  the 
newly  proposed  scheme  should  be  considered  at  the  level  of  the  fitting 
frequencies  appropriate  for  the  new  scheme.  The  frequency  dependence  of  the 
fluid  resistance  in  either  scheme  depends  on  the  typical  distance  between  the 
nodal  lines  of  the  imposed  displacements,  i.e.,  the  surface  expansion 
functions.  This  may  be  seen,  typically,  for  the  scheme  for  which  fitting 


results  are  shown  in  Figs.  1 to  7 . The  values  y (accession  to  mass)  were 

R 


fitted  for  one  frequency  above  and  one  frequency  below  the  peak  of  a curve. 
The  peak  shifts  to  higher  frequencies  as  the  mode  number  increases,  because 
the  nodal  distances  decrease  with  mode  number.  The  peak  is  located  at  quite 
a low  frequency  (f  ~ 0.3  kcps)  for  function  1,  higher  frequency  (f  ~ 1 kcps) 
for  function  5,  and  much  higher  frequency  for  function  12.  A shift  of 
fitting  frequencies  through  the  same  range  is  to  be  expected  for  the 
proposed  new  surface  expansion  functions.  The  range  f = 0.5  to  2.5  kcps 
was  therefore  studied  carefully. 


The  method  shown  in  the  Appendix  was  applied  to  determine  the  fluid 


resistance  (y_  and  y ) based  on  transformed  surface  expansion  functions  ip . ( s ) . 

R 1 .1 


Little  difference  in  the  importance  of  the  off-diagonal  terms  in  Y 


R 


(accession  to  mass)  was  observed  over  the  above  range  of  frequencies.  As 
typical.  Table  I shows  the  matrix  for  f = 1 kcps.  Inspection  of  the  Table 
indicates  that  the  off-diagonal  terms  are  nearly  always  less,  often  very 
much  less,  than  1/10  of  the  diagonal  ones.  The  ignoring  of  such  small 


terms  seems  reasonable. 
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The  check  on  the  importance  of  the  off-diagonal  terms  for  the  matrix  y 
(damping)  is  slightly  more  complex  because  of  the  nature  of  these  curves. 

They  have  the  character  of  Fig.  8.  The  values  are  practically  "ero  from  f = 0 
to  some  value  close  to  the  fitting  point,  (2.  As  a result,  the  prevailing  of 
the  diagonal  terms  is  of  importance  only  near  the  fitting  frequency  for  the 
respective  function.  At  frequencies  below  the  fitting  point,  the  diagonal  and 
the  off-diagonal  terms  are  several  orders  of  magnitude  smaller  than  near  the 
fitting  point.  The  situation  away  from  fitting  points  is  thus  of  no 
consequence.  This  is  illustrated  in  Tables  II  to  IV.  For  the  first  five 

A 

transformed  functions  ip.  the  fitting  points  are  below  or  up  to  f = 1 kcps,  and 

A 

Table  II  shows  the  matrix  y for  the  lowest  values  k,  j <5.  The  diagonal 

A 

terms  prevail,  not  so  strongly  as  for  y , but  still  with  factors  close  to  10. 
On  the  other  hand,  for  the  same  frequency  the  diagonal  terms  do  not  prevail 

for  higher  functions,  as  may  be  seen  from  Table  III  for  k,  j = 12  to  18. 

« -2  -3 

However,  the  values  y^  are  here  of  the  order  10  and  10  , i.e.,  extremely 

small  compared  to  the  diagonal  terms  in  Table  II,  which  are  of  the  order 

lO"*^.  For  values  k,  j > 12  the  fitting  frequencies  are  in  the  f = 3 to  4 kcps 

range,  and  Table  IV  shows  the  appropriate  portion  of  the  matrix  y . It  is 

seen  that  the  diagonal  terms  in  this  range  prevail  strongly,  with  factors 

larger  than  10. 

*) 

It  is  thus  concluded  that  the  proposal  to  disregard  off-diagonal  terms 
appears  to  be  a promising  approximation. 


It  is  noted  that  an  entirely  different  fitting  scheme,  involving  the  use 

of  convolution  integrals  was  suggested  at  a technical  meeting  at  DNA  in 

July  1974  by  Dr.  T.  Geers.  The  scheme  was  based  on  surface  expansion 

functions  which  lead  to  diagonalized  matrices  y and  y at  !2  = 0 and 

R I 

ft  ->  respectively.  This  scheme  also  ignored  off-diagonal  terms  at 
intermediate  frequencies.  The  degree  of  approximation  obtainable  by  any 
such  scheme,  shown  here  in  Tables  I to  IV,  has  however  not  been 
demonstrated  previously. 
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TABLE  I.  INERTIA  MATRIX  y BASED  ON  TRANSFORMED  SURFACE 
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FREQUENCY  IKCPS) 

FITTED  FLUID  RESPONSE  (REAL  PART) 
SURFACE  EXPANSION  FUNCTION  I,  N = 


□ fitting  point 


FREQUENCY  (KCPS) 


FITTED  FLUID  RESPONSE  (REAL  PART) 
SURFACE  EXPANSION  FUNCTION  2,  N = l 


FREQUENCY  (KCPS) 


FIG.  3 FITTED  FLUID  RESPONSE  (REAL  PART) 
SURFACE  EXPANSION  FUNCTION  5,  N=l 


FIG.  4 FITTED  FLUID  RESPONSE  (REAL  PART) 
SURFACE  EXPANSION  FUNCTION  12,  N = I 
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FREQUENCY  (KCPS) 


FITTED  FLUID  RESPONSE  (IMAGINARY  PART) 
SURFACE  EXPANSION  FUNCTION  I,  N=  I 


□ fitting  point 


FIG.  6 


FITTED  FLUID  RESPONSE  (IMAGINARY  PART) 
SURFACE  EXPANSION  FUNCTION  12,  N=l 


□ fitting  point 
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FIG.  7 FITTED  FLUID  RESPONSE  (IMAGINARY  PART) 
SURFACE  EXPANSION  FUNCTION  13,  N=l 


FREQUENCY 


FIG.  8 TYPICAL  BEHAVIOR  OF  A MAIN  DIAGONAL  TERM 

OF  DAMPING  MATRIX 
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APPENDIX 


Determination  of  Surface  Expansion  Functions  which  Diagonalize  the  Damping 
Matrix  for  High  Frequencies  and  the  Inertia  Matrix  for  Low  Frequencies. 


lb 


In  Ref.  [2],  elements  of  inertia  and  damping  matrices,  when  surface 

*) 


expansion  functions  are  used,  are  defined,  respectively,  as 


'kj 


P.(s)  ip  (s)  dA 
J K 


and 


= 


p.(s)  P.(s)  dA 
J k 


In  Eqs.  (A-l)  and  (A-2) , the  ^ are  the  surface  expansion  functions.  The 


replaced  by  the  , and  the  generalized  coordinates,  q . , of  Eq.  (123)  of 


Ref.  [2]  replace  the  q^  of  Eq.  (39)  of  Ref.  [2]. 


The  inertial  and  damping  components  of  the  generalized  forces 


corresponding  to  the  coordinates  q.,  defined  by  Eq . (126)  of  Ref.  [2],  are 


given  by  expressions  analogous  to  Eqs.  (110)  and  (98)  of  Ref.  [2]: 


Qlk  = l Ykj  qi 
.1 


52k  l Ykj  qj 

J 


(A-l) 


(A-2) 


p.  and  p^  are  defined  by  equations  analogous  to  Eqs.  (96)  and  (97)  of  Ref.  [2] 


In  Eq . (95)  of  Ref.  [2],  the  normal  components  of  the  shell  modes,  (p . , are 


(A- 3) 
(A-4) 


The  elements  of  the  inertia  and  damping  matrices  appearing  in  Eqs.  (A-l) 
and  (A-2)  are  related  to  the  elements  of  the  fluid  resistance  matrix  of 
Section  III  by  the  expressions 


Yk.i  - - nj  (Vki 


and  Ykj  = 


(YT).  . 
.1  I kj 


where  denotes  the  frequency  of  the  prescribed  time-harmonic  normal 


displacement.  Details  may  be  found  in  Ref.  [2] 
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If  the  normalizing  coefficients  p . of  Eq . (124)  of  Ref.  [2]  are  chosen 
equal  to  the  shell  surface  area,  the  orthogonality  condition  becomes 


ip.lp.  dA  = A 6.  . 
t J ij 


The  choice  of  an  orthogonal  set  of  surface  expansion  functions  is  not 
unique.  Given  one  set  ij;.(s),  new  functions  ip,  (s)  can  be  constructed  as 

1 K 

infinite  series  of  these: 


i|>k(s)  = l r^  ’ (s) 

i 


( A-  6 ) 


Because  the  problem  is  linear,  it  follows  from  their  definitions  that  the 

/s  £ 

pressures  p.(s)  and  p.(s)  corresponding  to  the  new  set  of  surface  expansion 

d A A 

functions  (J'j(s)  are  related  to  the  p^(s)  and  p^(s)  obtained  using  ip-(s)  by 
relations  similar  to  Eq.  (A-6), 

P,(s)  = l r.(i)  p.(s) 

J i 

p . (s)  = l r . ^ p . (s) 

i .1  i 


(A-8) 


Using  Eqs.  (A-5)  and  (A-6),  it  follows  that 


(i)  _ (j) 


ip.  dA  = A 7 r r 
l -j  £'  m m 

J m 


(A-9) 


Thus,  if  the  tj^  are  also  to  be  an  orthogonal  set  in  the  sense  of  Eq . (A-5), 


r (i)  r = 6.. 
mm  lj 


(A-10) 


As  was  shown  in  Section  J of  Ref.  [2],  the  use  of  any  set  of  orthogonal 
surface  expansion  functions  leads  to  a diagonal  damping  matrix  in  the  limit 
of  infinite  frequency,  i.e.. 
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11_fl  . ?W  * ' ‘ ' M ‘kj 

/V 

It  is  now  desired  to  choose  that  (unique)  orthogonal  set  of  ij/  [from  among 
all  those  obtainable  from  Eq.  (A-6)  with  coefficients  satisfying  Eq . (A-10)] 
which  diagonalizes  the  inertia  matrix  at  zero  frequency,  i.e.,  which 


(A-ll) 


results  in 


KJ  = 0 


b(k)  6, 


If  Eqs.  (A-6)  and  (A-7)  are  substituted  into  equations  analogous  to 

A A 

Eq.  (A-l) , with  y on  the  left-hand  side  instead  of  y , there  results 
kj  kj 

V = j F r (k)  r (j)  ? 
kj  Lt  n m nm 

m n 

Diagonalizat ion  at  zero  frequency  is  ensured  by  substituting  Eq . (A-13) 
into  Eq . (A-12) : 


U r (k>  r y ] 

n m nm  ~ ~ 

m n u - 0 


= b(k)  6. 


Thus,  the  set  (and  only  that  set)  of  surface  expansion  functions  which  is 
obtained  from  Eq.  (A-6)  with  expansion  coefficients  satisfying  both  Eqs. 
(A-10)  and  (A-14)  will  be  orthogonal  in  the  sense  of  Eq . (A-5)  and 
produce  diagonal  damping  and  inertia  matrices  at  infinite  and  zero 
frequencies,  respectively,  as  expressed  by  Eqs.  (A-ll)  and  (A-12). 


It  should  be  noted  that  the  elements  of  the  inertia  matrix,  y,  .,  are 

« kJ 

related  to  those  of  the  virtual  mass  matrix,  y . , by 

KJ 


n ->  o 


- \j  r“2  ■ Aj 


This  may  be  seen  by  comparing  the  definitions  of  Eqs.  (75)  and  (96) 
of  Ref.  [2]  with  surface  expansion  functions  used  instead  of  shell 
modes . 


(A-12) 


(A-13) 


( A- 1 4 ) 


20 


In  partial  matrix  form,  Eqs.  (A-6) , (A-10)  and  (A-14)  may  be  written  as 

5>.(s)  = J(s)  i(i)  (A-I5) 

r(l)  r(1)  = 6 . ( A- 1 6 ) 

- 3 J 

r(i)  j)_  r(j)  = b(i)  6 (A-17) 

fi  = 0 11 

Equations  (A-16)  and  (A-17)  constitute  a standard  eigenvalue  problem  in 
which  the  elements  of  the  (diagonal)  inertia  matrix  at  zero  frequency,  b^'\ 
are  the  eigenvalues  and  the  columns  of  expansion  coefficients  r^'^  are  the 
corresponding  eigenvectors.  The  solution  of  Eq.  (A-17)  which  satisfies 
Eq.  (A-16)  is  readily  seen  to  be 

0 - b(J)  I)  r(j)  = 0 (A- 18) 

j}=0 

in  which  I is  a unit  diagonal  matrix.  After  the  series  of  Eq . (A-6)  is 
truncated  by  considering  a finite  number  of  expansion  functions  tp . , the 
eigenvalue  problem  of  Eq . (A-18)  may  be  solved  by  any  standard  method. 


